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ABSTRACT 

Lens modeling of resolved image data has advanced rapidly over the past two 
decades. More recently pixel-based approaches, wherein the source is reconstructed on 
an irregular or adaptive grid, have become popular. Generally, the source reconstruc¬ 
tion takes place in a Bayesian framework and is guided by a set of sensible priors. We 
discuss the integration of a shapelets-based method into a Bayesian framework and 
quantify the required regularization. In such approaches, the source is reconstructed 
analytically, using a subset of a complete and orthonormal set of basis functions, 
known as shapelets. To calculate the flux in an image plane pixel, the pixel is split 
into two or more triangles (depending on the local magnification), and each shapelet 
basis function is integrated over the source plane. Source regularization (enforcement 
of priors on the source) can also be performed analytically. This approach greatly 
reduces the number of source parameters from the thousands to hundreds and re¬ 
sults in a posterior probability distribution that is much less noisy than pixel-based 
approaches. 
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1 INTRODUCTION 


Since the discovery of the first gravitationally lensed 
quasars, astronomers have used strong gravita¬ 
tional lensing to constrain cosmological parame- 
ters and explore distant quasars (s ee th e review by 
ISchneider. Kochanek fe Wambsganssl 120061 ) . Because 
quasars are typically unresolved, the positions, fluxes, and 
time delays of the observed images can be used to constrain 
the lens model. Galaxies, on the other hand, are extended 
objects, and there is quite a lot of information buried in 
how the surface brightness of the images varies pixel to 
pixel. 


Early attempts to model the source’s (the galaxy’s) 
light profile assumed elliptica l symmetry and analyzed 
isoph otal shapes (e.g., iBIandford. Surpi fe Kundid 
200ll) or peak surface brightness curves (e.g., 


Kochanek. Keeton fe McLeodl l200ll ) in Einstein rings. 
More general strategies used an iterative, two-loop 
method to solve for the best lens and source model s 
llWallington, Kochanek fe Narav^ ll99^jKoqpmnndl20d^). 


As this can be time-consuming, Warren fc D ye (2003) 
introduced a semilinear process in which, for a given lens 
model, the source could be solved for analytically, greatly 
speeding up the process. 


More recent methods, however, reconstruct the source 
on an irregular, p i xelate d grid dDve fe Warrenl 120051 : 
IVegetti fe Koopmand 120091 ). Doing so places tighter con¬ 
straints on lens models and can allow substructure 
to be detected as well ( e .g., 


_ Suvu fe Halkola 

Ve getti. Czoske fe Koopmansl 1 201C ; IVegetti et al. 


201C : 


2010 . 


20121) . Additionally, there has been an increased interest 


in exploring the structure of the galaxies being lensed, 
and it is possible to recover small-scale features in the 
de-lensed source ( e.g., ISharon et al.ll20lT iDve et ahll2013l : 
iRvbak et al.ll2015l) . 


However, pixel-based reconstruction algorithms require 
a large number of source parameters, which can reach into 
the thousands. Even on modern computers, it can take 
weeks to months to analyze a single gravitat i onally lensed 
system. Recently, iBirrer, Amara fc Refregien (120151 ) intro¬ 
duced a shapelets-based method, wherein the source is re¬ 
constructed analytically. Shapelets are a family of functions 
that are mathematically complete and orthonormal, and the 
source’s surface bright ness is expanded using a a finite sub¬ 
set of the shapelets. As lRefregierl d2003l ) have demonstrated, 
the Hermite polynomials can be used to construct shapelets 
that can reproduce galaxy shapes using a small number of 
shapelets. In practice, only a few hundred source parameters 
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are needed to account for the surface brightness distribution 
seen in the observed images. 

Here, we extend the analysis by formulating the 
shapelets-based source reconstruction in a Bayesian frame¬ 
work. In order to stabilize the source reconstruction, we 
place priors on the source’s surface brightness (i.e., regular¬ 
izing the source) and optimize the strength of the regulariza¬ 
tion by maximizing the Bayesian evidence. Additionally, we 
present an image-plane tiling algorithm and integrate the 
shapelets over the source plane, as would occur naturally 
during an observation. We then compare results from the 
shapelets-based method to existing pixel-based methods. 


2 SHAPELETS 

2.1 One-dimensional basis functions 

In order to describe a galaxy usin g shapele ts , we adopt 
the basis functions deve loped in iRefregieu (120031 1 and 
iMassev fe Refregied d2005ll . Here, only the required back¬ 
ground information to keep this work self-contained is given. 
For a more complete explanation, please see the aforemen¬ 
tioned pa pers. _ _ 

From [Refregied ( |2003l ) the one-dimensional, dimension¬ 
less, Cartesian basis functions are given by 

bn{x) = (2 n 7r2n!) 5 H n (x) e 2", (1) 

where H n is the n th order Hermite polynomial. For describ¬ 
ing galaxy light profiles, the basis functions are scaled to 
obtain 

(j> n (x) = (2) 

where /3 is the characteristic length scale of the galaxy. Both 
the dimensional and dimensionless basis functions are, sepa¬ 
rately, complete and orthonormal. Thus a galaxy, no matter 
how complex, can be described using these shapelets. 


IMassev fc Refregied (l2005l l the two-dimensional basis func¬ 
tions in this case are given by 


Xn.,m(x, 9 ) — 


(- 1 ) 


Hml 


( «Hro| )| x i 




r |m| x 


(6) 


- M 


e e 


where r and 9 are the radial and angular coordinates, re¬ 
spectively, f3 is (as before) the characteristic scale length, 
and Lp(x) are the associated Laguerre polynomials. The 
order of the basis function is determined by n and m. n 
can be any non-negative integer, m can take on values be¬ 
tween —n and n but must increase in steps of two. That is, 
m = —n, —n + 2,..., n — 2, n. 

As before, the polar shapelets are complete and or¬ 
thonormal, and we can expand the source’s surface bright¬ 
ness as 


S(r,9) = Y^ c n,™Xn,m(r,0;P), (7) 

71=0 m=—n 


and the shapelet coefficients are given by the overlap integral 

Cn,m = j j rdr d9 S(r, 9) Xn,,m(r, 9; /3) (8) 

As mentioned, the summation over m assumes m increases 
in steps of two. 

Unlike the Cartesian case, the polar basis functions and 
their shapelet coefficients are complex. The magnitude of a 
given complex shapelet coefficient determines how strongly 
that shapelet contributes to the surface brightness; the phase 
determines the angle of orientation. We note that although 
the complex polar coefficients contain twice as many free pa¬ 
rameters than do an equivalent number of Cartesian coeffi¬ 
cients, the polar coefficients are not all independent. Specif¬ 
ically, Cn,m = Cn,-m, where the asterisk denotes complex 
conjugation. This relation implies that there is the same 
amount of information encoded in both the Cartesian and 
polar coefficients. 


2.2 Two-dimensional, Cartesian shapelets 

Two-dimensional, Cartesian basis functions can be con¬ 
structed from the tensor product of the one-dimensional 
functions. We have, 

^n x ,n v (x,y) = (3) 

For a galaxy, its surface brightness distribution, S(x , y), 

can be expanded as 

oo oo 

S(x,y)=^2 ^2 c ^^ y ^n X: n y (x,y), (4) 

n x =0 riy =0 

and the shapelet coefficients, c nx , ny , are given by the overlap 
integral 

dxdyS{x,y) $„ x , ny {x, y). (5) 


2.3 Two-dimensional, polar shapelets 

We can also describe the source using a polar coordinate 
system, which is a more natural choice for galaxies. From 


2.4 Practical issues 

There are several free parameters when perfor ming a 
shapelets reconstruction. IMassev fe Refregien (120051 1 present 
a set of criteria and an iterative approach for optimizing the 
shapelet center and scale, as well as the maximum order of 
the shapelets basis functions. Because this approach is time- 
consuming when searching the lens model parameter space, 
we instead fix all three of these parameters. The maximum 
order of t he shapelets is fixed but can be adjusted later if 
necessary. iBirrer, Amara 8z Refregierl (|2015l l also make these 
simplifying assumptions and apply their method to the grav¬ 
itational lens RXJ1131-1231. In H3.2I and ( 13.31 we discuss 
issues related to integrating the shapelets over the source 
plane and fixing the maximum shapelets order. We also ex¬ 
plore the effects the shapelets order can have on the recon¬ 
struction in 0 

The shapelets center and scale, on the other hand, vary 
from one lens model to the next. To estimate the center, 
we ray-trace image pixels with surface brightness values 
greater than three times the noise level. The surface bright¬ 
ness weighted positions of the subset of image pixels are used 
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to fix the shapelets center. Similarly, the surface brightness 
weighted distances of the ray-traced pixels to the center are 
used to fix the scaleQ 

In order to easily integrate over the the shapelets basis 
functions (see H3.2I) , we choose to use the Cartesian shapelets 
when performing the reco nstruction. However, as noted in 
iMassev fc Refregierl (20051 ). these can be later converted to 
polar shapelets, which are useful for obtaining analytic func¬ 
tions for properties of the source, such as the azimuthally 
averaged radial light profile. 


be the sum of the squares of the derivatives of the source’s 
surface brightness, within a factor of two. 

Minimizing P leads to an equation for the reconstructed 
source 

s r = F -1 D, (13) 

where F = L T CJ 1 L + AH T H and D = L T Cj 1 d. 

We note that it is also possible to non-linearly solve for 
the optimal regularization strength by minimi zing P. For 
detail s of the derivation and optimization, see ISuvu et al.l 
< 20061 ) . 


3 SOURCE RECONSTRUCTION 
METHODOLOGY 

With a few modifications, shapelets can enter into the estab¬ 
lished Bayesian frameworks for analyzing resolved gravita¬ 
tionally lensed images using pixelated source grids. In these 
frameworks, the ideal, reconstructed source reproduces the 
data, while not violating any priors placed on it. In this 
section we first briefly review the formal Bayesian frame¬ 
work, which has been discussed in detail by Warren fe Dvd 
(20031).I Suvu et ahl (l2006| ). IVegetti fe Koopmansl ( 20091 ) . and 
iBrewer fe Lewisl ~l 20061 ). Then, we show how shapelets can 
replace a pixelated source grid. We note that lens model 
ranking is not discussed here but is addressed in the afore¬ 
mentioned papers. 


3.1 Bayesian framework 

In the absence of scattering or absorption of light, lensing 
conserves surface brightness. We can therefore relate the 
source and lensed images as 

d = Ls + n, (9) 

where L is a linear “lensing operator” that include grav¬ 
itational, atmospheric, and instrumental effects, s and d 
are vectors containing the surface brightness values in the 
source plane and image plane, respectively, and n is the noise 
present in the data. 

To reconstruct the source, a penalty function, P, is in¬ 
troduced: 


P = G + \R, (10) 

where G quantifies the mismatch between the data and the 
model, R quantifies how much the source violates priors 
placed upon it, and A is the regularization strength, which 
determines how strongly the priors are enforced. Specifically, 

G=i(Ls-d) T CJ 1 (Ls-d), (11) 

where Cd is the noise covariance matrix of the data, and 

R=i(Hs) T Hs, (12) 

where H encodes the priors placed on the source. For exam¬ 
ple, H could act on a source vector s to produce a vector of 
derivatives of the source’s surface brightness. R would then 


3.2 Constructing the lensing operator 

The primary difference between pixel-based and shapelets- 
based source reconstruction methods lies in the construction 
of the lensing operat or. In pixel-based methods (see, e.g., 
IVegetti fe Koopmandl2009l ). each image-plane pixel is ray- 
traced to the source plane. The irregular, pixelated source 
grid is triangulated, and the three source pixels that enclose 
the ray-traced data pixel are identified. Bilinear weights are 
placed on these three source pixels so that the total surface 
brightness matches the data pixel’s surface brightness. In 
this way, the lensing operator is constructed, and there are 
at most three entries per row0 

When using shapelets, there is no source grid, and the 
source vector does not contain surface brightness values. In¬ 
stead, the source vector contains the strengths of the basis 
functions, and the position along the source vector deter¬ 
mines the order of the basis function. In order to create 
the lensing operator, each data pixel is ray-traced to the 
source plane, as before. The location of a ray-traced pixel de¬ 
termines how strongly each basis function contributes (i.e., 
the basis function evaluated at the location) to the surface 
brightness of that particular data pixel. Thus, the lensing 
operator contains no empty elements. However, this does 
not necessarily imply that the shapelets-based approach is 
slower or more computationally expensive, as the dimen¬ 
sions of the lensing operator are much fewer than in the 
pixel-based case. 

In the absence of a gravitational lens, there is a mini¬ 
mum (set by the pixel size) and maximum (set by the image 
dimensions) sca le over which information can be obtained 
(Refregierl 120031 ). These extremes will constrain the maxi¬ 
mum order of the basis functions and shapelets scale that 
should be used in reconstructing the image. However, be¬ 
cause gravitational lensing magnifies some regions of the im¬ 
age more than others, determining the maximum shapelet 
order and scale is not straightforward. Whereas image pix¬ 
els in highly magnified regions may require higher-order 
shapelets, pixels in lower magnification regions may, essen¬ 
tially, randomly sample the highly oscillatory shapelet. 

In order to combat this issue, when calculating the sur¬ 
face brightness in an image pixel, we integrate each shapelet 
over the ray-traced region in the source plane. Because 
highly oscillatory shapelets will tend to average to zero over 
large regions, we fix the maximum order of the shapelets at 
some reasonable value (i.e., 20 x 20 — 50 x 50). The shapelet 


1 If the noise is not uniform across the image, inverse noise 2 Before convolution with the instrumental PSF or any other 

weighting can be used as well. kernels. 
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integration is complicated by the fact that a square image 
pixel will not map to a square region in the source plane, 
and it may, in fact, overlap on itself. We thus split each im¬ 
age pixel into N sp nt right triangles, because triangles will 
always map to triangles and are more easily handled 0 For 
a particular image plane pixel, each of the JV sp iit integrals 
will give the average surface brightness over that triangle. 
The average surface brightness over the image plane pixel p 
is then given by 

^split „ .. 

s P = N-\ it E ^ _1 EE J j dxdy$n X: n y (x,y), (14) 

t—1 77. X 71 y 

where At is the area of the ray-traced triangle. 

IVspiit can vary from one image pixel to the next. We 
choose the value based on the local magnification, so that, 
for a particular pixel, 


N spm =l2p\, (15) 

where p is the magnification at the center of the pixel and 
the half square brackets denote the floor function. This pre¬ 
scription allows iVspiit to follow the magnification relatively 
smoothly across the image plane. However, we allow for min¬ 
imum (two or greater) and maximum values of AT sp i it to be 
enforced. 

Additionally, a ray-traced triangle will not remain a 
right triangle, and the integration will be performed in two 
steps as follows. We choose to fix the lower and upper lim¬ 
its of integration in the x direction. Subsequently, the lower 
and upper limits of integration in the y direction will be 
a function of the x-coordinate. We therefore split each ray- 
traced triangle into two smaller triangles, such that both tri¬ 
angles share the vertex with the median ^-coordinate, and 
the line segment that splits the triangle has infinite slope. A 
schematic diagram is shown in Fig[T| 

There is no closed-form solution for integrating the tri¬ 
angles in the source plane. We instead present a series of 
recursion relation which can be used to quickly calculate the 
integ rals. For fixed limits of integration. iMassev fe Refregien 
(120051) show that the one-dimensional integral over the 
shapelets is given by 



Figure 1. Schematic diagram illustrating the shapelets integra¬ 
tion. First, the image pixels are split into right triangles, where 
the number of splittings depends on the local magnification (or 
can be fixed to a constant value). These triangles are then mapped 
back to the source plane. Focusing on one ray-traced triangle, it 
is further divided into two smaller triangles, such that the two 
share the vertex with the median ^-coordinate, and the line seg¬ 
ment that splits the triangle has infinite slope. Then, two integra¬ 
tions are performed. Each integration has fixed upper and lower 
limits over the x-coordinate. The limits of the integral over the 
^/-coordinate are, however, functions of the ^-coordinate and are 
denoted by the solid, blue and red lines. 


by 


d rri2X-\-b2 


l n x ,riy — 


/ dx I d y <bn x (x) 4> ny (y) 

c m\x-\-bi 
d 

J dx <j>n x (x) I„ y (mix + bi,m,2X + b 2 ) 


n y 1 r it 

‘ n x ,riy —2 i Jn x ,n y — 1 1 


where we have defined 


0 

I„(a, b) = J dx'c/)n(x') 

a 

= j 71 1 In. — 2 (u, b ) - P\f — <j) n -l(x') 

V n V n 


(16) 


We seek to evaluate the two-dimensional integral given 


3 This i mage-plane til ing scheme is analogous to the one pre¬ 
sented in lKeetonl (1200ll) , although differences are discussed in the 
text. 


’’lx i n y 



rri2X+b2 
x' =m\x+bi 


(18) 

Note that x' is evaluated as a function of the dummy variable 
x over which the integral is evaluated. 

Before calculating Jn x ,n y , we first present recursion re¬ 
lations for the In x ,n y with n y = 0,1. From hereon out in 
order to keep the equations tractable, we will no longer ex¬ 
plicitly write the evaluation of dummy variables at the limits 
of the integrals. For example, we will denote an expression 
of the form (mx + 6)|™=m i; b=b 1 l*= c = ((m 2 x + b 2 ) - {rmx + 
bi))\i- c = ((m 2 d+fo 2 ) —(mid+&i)) —((m 2 C+fc 2 ) —(rmc+fei)) 
simply as (mx + b). Thus, the remainder of the equations 
in this section may actually consist of up to four times as 
many terms as are explicitly written. Using this notation, 
we have 
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above relation and the following: 


and 


(19) 


(m 2 + l)i/h^ 

-1 




Jo ,n y — 

(m* + l)y/n v + 1 


In x ,1 — 

rinl r ti- 1 i / \ -(mx-\-b) 2 / (28 2 ) 

2 / 32 [tt 4 ] 0 n:B _i(x)e v ^ ; /K p ’ 
/3 

- (m 2 - 1 )Vn x - ll nx - 2 ,i 


— 4>n y -i(mx + b)e x 7(2,3 ) 
7T4 ./nZ 


+ %/ 26 / /3Jo,n y — l 

+ (m 2 — 1 )(n y — l)/^ThJJo,n y -2 


(25) 


fix tn X: o — tt 


I (mx + b\ , . . 

p2eri {-^vrr n *- i{x) 


- ml nx - 1,1 

+ \/nx — l/n x _ 2 ,o 


( 20 ) 


r _ V n v + 2 r 

— ^O.n-u+l 

m y 

n y T V2b T 

+ m /TTT °’" rl o J0,riy , 

Hiy/ I by _L \llfj 


(26) 


The starting points for these relations are given by 

= P x + m{b + mx) \ ^- 6 3 /(2 g2 (1+m 2 )) 

%/l + m 2 \ a/ 2/3 2 (1 + m 2 ) / 


ii.o = — /Serf 


ft + rnx^ -x 2 /{ 2 p 2 ) 


( 21 ) 


^2,5 


+ 


_£m_ prf Y x + m(b + mx) ^ ^-b 2 /(2/3 2 (i+m 2 )) 
%/l + w- 2 V 2/3 2 (1 + m 2 ) / 


( 22 ) 


and 


7l r =e -(i> 2 +26m a; +(l+m 2 )x 2 )/(2/3 2 )^ 1 + m 2s-3/2 ^1/2 x 


2/3 \/1 + m 2 + 72^6me (:c+m(6+mx))2/(2/32(1+m2)) x 


erf 


x + m(6 + mx) \ 
V2/9 2 (l + m 2 )A 


(23) 

To our knowledge, there does not exist a closed form for 
/o,o- It can be computed numerically, or by approximating 
the error function as the product of an exponential function 
and a truncated power series expansion, 7o,o can be com¬ 
puted analytically to high accuracy. See the appendix for 
details. 

It remains to evaluate the integrals Jn x ,n y ■ Again, we 
present recursion relations to compute these quickly. 


J n x ,riy 


(m 2 + 1 )y/nH 


2/3 2 


p <j>n x -i{x)<j> n (mx + b) 


\J Tly T 1 

- bmp -1 V2J nx -\, ny 

(tu 1 'jy/fix fJ r n x —2,n 
2 mn v 


Jn x —l,n v — l 


(24) 


+ 1 

All the Jn x ,n y integrals are fully specified with the 


Jn x ,0 — In x ,l , 


(27) 


and 


do,i = 


Z i 1 

772 +1 


m/3e 


-((mx+b) 2 +x 2 )/(2/3 2 ) 


- ( m 2 + l)-he- b2 ^ 2+1 ^eri( ^ ± 1)x ± mb 

V \j2(m? + l)/3 


(28) 


3.3 Regularizing the source 

Just as with pixel-based source reconstruction methods, 
regularization is necessary, especially since the maximum 
shapelet order is fixed and not set by the data. Several 
forms of regula rization have bec ome popular in the liter¬ 
ature (see, e.g.. ISuvu et al.ll2006h . We present formulae for 
calculating the zeroeth, first, and second derivatives of the 
surface brightess and for calculating the RMS size of the 
source. Shapelet indices with a hat above them (e.g., fix) 
denote the maximum order of the shapelet along the speci¬ 
fied coordinate. 

Zeroth order regularization penalizes sources with large 
surface brightness values and is given by the integral over 
all space of the square of the surface brightness: 

oo oo 2 

r o = - J2 j J dxd+Cn<Mx)j ( 29 ) 

n ' ' n 


If we instead want to minimize the change in surface 
brightness, then a possible choice for gradient regularization 
is given by the integral over all space of the square of the 


















































6 Tagore & Jackson 


magnitude of the gradient of the surface brightness: 




did yl c n V<f> n (x) 


— OO —OO 


n x + l 


4 « 2 'E ( Cl ’™y) + E (gfnjcn*- l,n„) 

77 y 77 x —n x 

fix —1 

“1“ ^ ^ ^y/TlxCn x — l,ny 'V ^x ^- c n x + l,n y ^ 

n x = l 

^ 2 ^y+1 

+ 4«2 £ ( C "x,l) + E (\/ny C n,x,riy-2) 

Tl X Tly'^fly 

fly-\ 

+ 'y \ (yy/TlyCnz.,riy — 1 n y T l c n X! n y + l^ 


(30) 


If we want to ensure that the surface brightness is 
smoother over even larger scales, then a possible choice for 
curvature regularization is given by the integral over all 
space of the square of the Laplacian of the surface bright¬ 
ness: 


R2 = ~^2 / didy( c n V 2 <f>n(x) 


1 


— OO —OO 

1 r 1 


E E {K lx (n x ,n y ) + K ly (n x ,n y )) 2 


n x =0 L n v =0 


8/3 4 

77 ry 2 

H - ^ ^ (-^1 x(jlx i Hy') H" -^2y {jlx i ^7/)) 
n y =2 

fly 

H - ^ ^ (-^1 xijTlx •> Hy') H - -^3y {j^x -> ^t/)) 

n y =fi y — l 

n x 2 r- 1 

+ EE {^K2x{jl j x Tly') “I - -^ly (jT x ? Tly)) 

Tljj —2 77 y — 0 

rty 2 

H - ^ ^ (-^2 xiyixi'Hy') + ^2y {j^x ? Tly)) 

77y-2 

fly 

+ ^ (K 2x (n x ,n y ) + K 3 y(n x ,n y )) 2 

Tly -77 y 1 

r 1 

+ E E {K 3x (n x ,n y ) + K ly (n x ,n y )) 2 

rh x — TL X 1 ?7y —0 

rty 2 

+ ^2 (K 3x (n x ,n y ) + K 2y (n x ,n y )) 2 

TLy -2 

fly 

+ ^ (K 3x (n x ,n y ) + K 3 y(n x ,n y )) 2 

Tly -77 y 1 

n x +2 _ 

+ E E (\/n x {n x - l)c nx - 2 ,n.y) 2 

n y =0 n x =fi x +1 
77;c nj,+2 

+ E E ( Vn y (n y - 1) Cn x ,n y —2) 

77 3 ; =0 Tly = fiy + 1 


(31) 


where 

^la:(^a:) %) = ijlx H“ l)(^a; “I - 2)"c n;:c _|_2,77 y (2?7.^ H - l)Cn cc ,77y 
K ly (n x ,n y ) = v / (ETl)(n7+2)cn a ,,^ + 2 - ( 2 n„ + l)c na; ,„ w 
K 2 x (n x ,n y ) = ^ (^ + l)(nx + 2 )c„ x+2 ,„ y 

(2?7.x “1“ l^n^ny “f \/Tlx^Tlx l)Cn x —2,n y 

K 2y (n x ,n y ) = -y/ ( n y + l)(n,, + 2)c„ Xiny + 2 

— (2 Uy + + 'J'n-yijly — V)Cn Xt n y — 2 

K 3 x(n x ,Tiy') — (2Ux “{“ l)Cn x ,7i y "F \/Tlx{ri X l)Cn x —2,n y 

K 3y (n x , n y ) = —(2 n y + l)^^™,, + \Zn y {n y — l)cn x ,n y —2 

(32) 

Finally, if the source is known or expected to be com¬ 
pact, we can instead minimize a sort of un-normalized RMS 
radius, given by 


OO OO 

Rsize = 7T E J j dxdy^y/x 2 + y 2 C n $n(x)^ 

n — 00—00 

T/ , 2 

•E («•-.) + E 


r 

4 


/ Tlx Cria; — 1 


>riy) 


+ ^ ^ ^\JTlx Cn x — l,riy H - \/^x +l,77y ^ 

77 a; = 1 

2 T 2 77y + l 

+ — E/ ( C "a:. 1 ) + E/ 

fix L n y =fi y 

™y 1 _ 2 

+ (^\/nyCn x ,ny -1 + \/riy + lCn x ,n y +2j 


(33) 


4 TESTING THE SHAPELETS 

4.1 Varying shapelet parameters 

Here, we look for effects of varying the maximum shapelet 
order and the maximum value for V S piit on the source recon¬ 
struction. Using mock data exhibiting complex morpholo¬ 
gies, we also test the ability of the shapelets reconstruction 
to capture both small and large scale features. The source 
being lensed is comprised of two merging galaxies, with a 
bridge connecting the two, and crosses a caustic, so that 
some regions are doubly imaged, while others are quadruply 
imaged. The lens is a singular isothermal ellipsoid with ellip- 
ticity of 0.3 and position angle of 60°, east of north (EoN). 
The Einstein radius along the major axis is 3”, where the 
pixel scale is 0.05 arcs/pixel. An elliptical Gaussian PSF has 
been applied, with a full-width half-maximum (FWHM) of 
O^SO” and 0.165” along the major and minor axes, respec¬ 
tively, and a position angle of —55° EoN. Finally, additive 
Gaussian noise has been applied, so that the final data has 
a peak S/N« 21. With the lens model fixed at the true 
parameters, we perform source reconstructions for several 
different shapelets configurations, shown in Fig. 0 We find 
that using 10 x 10 basis functions is not adequate to cap¬ 
ture all of features in the data. 20 x 20 basis functions or 
more do a reasonable job of capturing all features, with the 
50 x 50 case resulting in a reduced y 2 = 1 (y 2 per pixel). For 
this particular mock data, increasing the maximum value of 
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Figure 2. Shapelets reconstruction for several shapelets configurations. Top row, left to right: data, source reconstruction, true source, 
model residuals. 50 X 50 shapelets basis function were used. Bottom row, left to right: Model residuals for runs with 10 X 10, 20 X 20, 
30 X 30, 40 X 40 shapelets basis functions. All runs used a maximum ,V sp | P of 2 X 10 6 and used optimal gradient regularization. We find 
that, for this particular mock data, lowering A r sp i; t has no significant effect. 


lV S piit from two to 2 x 10 6 does not result in a significant 
improvement for any of the model residuals shown. 

However, we find that N sp lit can become important if 
the source size is reduced and/or moved near a caustic. In 
order to explore this further, using the same lens as was 
used in the preceding mock data analysis, we fix the source 
to be a circular Gaussian. The source position is fixed 0.25" 
along the major axis (within the quadruply imaged region), 
while the FWHM of the source and A r sp iit are varied. We 
then limit the shapelets to 1 x 1 basis functions (identical to 
a circular Gaussian) and fix the shapelets center and scale 
to the true values for each run. Using this setup, we can 
attribute any inability of the shapelets integration to recover 
the true image plane flux to an insufficiently small value of 

Wpii t.. 

We find that, for larger source sizes, IV S piit has little to 
no effect, but as the source size becomes smaller, a higher 
JVgpiit is needed to accurately calculate the flux. See Fig. [3] 
For the smallest source size (3.5 x 10 -3 pixels at FWHM) 
and highest value of iV S piit (8.4 x 10 6 ), the image plane flux 
of each image calculated using the shapelets method agrees 
with the results for a point source to within 2%Q For the 
next smallest source, (3.5 x 10~ 2 pixels), agreement is to 
within 0.5%. For larger sources, finite source effects become 
important. For these reasons, we quote the percent error for 
each source size relative to the flux calculated using the max¬ 
imum value of IVspiit, which is limited by the computational 
resources required, for that particular source size. 

For radio observations of quasars, finite source effects of 
the emitting region can become significant and distinguish¬ 
able from a point source model. In such cases, Wplit can 
play an important role in accurately recovering image-plane 
fluxes. In Jackson, et. al. (2015, in prep.), we investigate 
such finite source effects from several radio-quiet quasars. 


4 As the Gaussian source become smaller, the total flux calcu¬ 
lated should approach the analytic results for a point source. 



Figure 3. Fractional error in the flux calculation for several dif¬ 
ferent source sizes as a function of N sp n t . Different lines represent 
different source sizes, measured in pixel units at FWHM. The frac¬ 
tional error is computed relative to the flux calculated using the 
maximum A r sp | Pi for that source size. Note: a fractional error of 
-1 indicates that no flux was recovered in the image plane. 

4.2 Varying lens model parameters 

Here, we turn to examining the ability of the shapelets 
method to infer lens model parameters and compare the re¬ 
sults against the well-established pixel-based approach. The 
lens model is the same as in the previous section, but the 
source is now an elliptical Gaussian, placed in a cusp configu¬ 
ration, and additive Gaussian noise is applied, so that a peak 
S/N~ 41 is achieved. In Fig. [4] we show a one-dimensional 








































8 Tagore & Jackson 




ellipticity 

Figure 4. One dimensional projection of — 21n£. All model pa¬ 
rameters are fixed at the true values, while ellipticity is varied 
about its true value of e = 0.3. Top: Projection using shapelets, 
with 20 X 20 shapelets basis functions and a maximum Al sp iit of 
~ 3 X 10 4 . Bottom: Projection using fully adaptive grid, twice the 
number of image pixels as are source pixels. Both panels use cur¬ 
vature regularization, with the regularization strength optimized 
at each step. 

projection of the evidence as a function of ellipticity. All 
other lens model parameters were fixed at their true values. 
Because shapelets do not require a source grid and because 
the pixelization of the image plane is accounted for by in¬ 
tegrating the shapelets in the source plane, the resulting 
X curve is relatively smooth. We find typical fluctuations 
of Ay 2 < 1 using the shapelets-based method, while pixel- 
based methods can result in Ay 2 ~ 10. 

Finally, we also perform a full lens model parameter 
space exploration, varying the Einstein radius, position, el¬ 
lipticity, and position angle of the lens. The source is still an 
elliptical Gaussian in a cusp configuration but, this time, the 
data have a peak S/N« 300 and an elliptical Gaussian PSF 
has been applied, with FWHMs of 0.236 ,/ and 0.165 ,/ along 
the major and minor axes, respectively, and a position an¬ 
gle of —55° EoN. Joint posterior probability distributions for 



(Einstein radius - 3) x 10 4 [arcs] 



Figure 5. Joint posterior probability distributions for (top) Ein¬ 
stein radius and ellipticity and (bottom) right ascension and dec¬ 
lination. For both panels, the true values of the lens model pa¬ 
rameters have been subtracted and scaled, so that the origins of 
the coordinate systems coincide with the true lens model param¬ 
eters. Because of the high peak S/Nra 300, the uncertainties on 
the parameters are small, but (not unexpectedly) all source re¬ 
construction methods show a clear bias on the estimate for all 
parameters, with the exception of the Einstein radius. 

two pairs of lens model parameters are shown in Fig. [5] The 
shapelets-based method does a reasonable job of inferring 
the lens model parameters. Except for the right ascension, 
all parameters are recovered at the 95% confidence level. The 
Einstein radius is the most robust of all the parameters, for 
all source reconstructions methods. However, there is a clear 
bias in the positional parameters and the ellipticity. For data 
of high S/N, su ch as these, an observed bias is not unex¬ 
pecte d (see, e. g, iTagore fe Keetoij|20l4 iNightingale fe Dvcj 
120141) . Interestingly, the biases in the shapelets and pixel- 
based approaches seem to be in opposite directions, but we 
have not explored this result fully. 


5 CONCLUSIONS 

We have developed an implementation of a shapelets- 
based, as opposed to pixel-based, approach for reconstruct- 
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ing resolved, gravitationally lensed images. Expanding the 
source’s surface brightness as a finite series of shapelets 
requires a choice of several parameters. We have cho¬ 
sen to fix the center and scale of the shapelets using 
the data and lens mode l them selves. Using mock data, 
iBirrer. Amara fc Refregief feoisl ') are able to accurately re¬ 
construct their input source model and fit the lensed im¬ 
ages. We find that the source reconstruction can, for many 
lens configurations, lead to unphysical source reconstruc¬ 
tions. We therefore fix the maximum shapelets order and 
the maximum number of image pixel splittings N sp \it before¬ 
hand and regularize the source using a Bayesian framework, 
which stabilizes the reconstruction. If there are not enough 
shapelet coefficients to accurately describe the source, model 
residuals will clearly indicate this. We also find that for more 
extended sources, A^piit does not play a significant role. For 
smaller sources, Af sp iit can be increased to test for robustness 
of the source reconstruction. 

The advantages of the shapelets lie in the relatively 
small number of parameters needed to reproduce the data. 
We find that even for complicated morphologies, 20 x 20 = 
400 source parameters are sufficient to capture most details. 
On the other hand, creating the lensing operator requires 
evaluating every shapelet coefficient for every pixel, which 
results in a dense, completely filled lensing operator matrix. 
In practice, this does not result in a decrease in performance, 
because the recursion relations presented in 1 )3.21 allow for 
quick calculations as long as jV ap ii t is not unnecessarily large. 
Moreover, a bottleneck in pixel-based algorithms occurs dur¬ 
ing optimization of the regularization strength[f] Because of 
the small number of source parameters, this step can be 
performed quickly, giving an unbiased estimator of the lens 
model being ranked. 

We note, however, that given data exhibiting more com¬ 
plicated features than we have explored here (e.g., many 
starburst regions that are individually resolved or multi¬ 
ple, widely separated source galaxies), a large number of 
shapelets may be required, resulting in a decrease in per¬ 
formance. A possible solution to this is to create multiple 
shapelet expansions centered at the appropriate positions, 
but we have not explored this possibility. 


APPENDIX A: INTEGRAL APPROXIMATION 

The 7o,o integral cannot be calculated analytically. It can be 
reduced to a single integral: 

c 

(md+b)/(V2p) 

- t f duerf(u)e- (v ^“- b)2/(2/32m2) . 
m J 

(mc+i>)/(V2/J) 

By approximating the error function as the product of an 
exponential and a truncated power series expansion, the in- 


5 Some authors omit this optimization and account for its effects 
in other ways. 


tegral can be computed analytically. Specifically, we use 
eri(u) l — e dlu+d2U (1 + d 3 u + d^u 2 + d$u 3 + deid). (A2) 
Optimizing the six free parameters, we find that 
di = -1.6093661701173445, 
d 2 = -0.9144569742603699, 
d 3 = 0.48105948656562003, 
d 4 = 0.3927744807707903, (A3) 

d 5 = 0.05250001720254402, 
and 

d 6 = 0.036174940195734834. 

Although the approximation is valid only for positive 
arguments, negative arguments can be evaluated by real¬ 
izing the error function is odd. Over the entire domain of 
the function, we find a maximum absolute error of 2 x 10~ 6 
and maximum relative error of 6 x 10 -5 . Plugging this ap¬ 
proximation into eq. IATI an analytic solution can be found, 
although because of the large number of terms in the result¬ 
ing expression, it cannot be sensibly reproduced here. 
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